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1.  INTRODUCTION 

This  report  covers  progress  made  over  the  entire  nine  month  period  of  the  Phase  1  STTR  project  en¬ 
titled  High- Order  Modeling  of  Applied  Multi-Physies  Phenomena^  awarded  to  Scientific  Simulations 
LLC  and  performed  in  collaboration  with  the  University  of  Wyoming. 

The  objective  of  this  project  has  been  to  demonstrate  the  potential  of  a  novel  simulation 
capability  based  on  high-order  discretizations  in  both  space  and  time  for  application  to  practical 
engineering  problems  involving  complex  physical  phenomena  and  complicated  geometries.  The 
long-term  goal  is  to  develop  a  tool  which  can  accurately  handle  multiphysics  simulations,  both  in 
analysis  mode,  and  for  design  optimization  purposes. 

The  general  approach  used  relies  on  high-order  (up  to  6th  order)  Discontinuous  Galerkin  dis¬ 
cretizations  in  space,  and  backwards  difference  as  well  as  higher-order  implicit  Runge-Kutta  tem¬ 
poral  discretizations  (up  to  5th  order).  The  work  performed  in  this  project  is  intended  to  show  how 
the  favorable  asymptotic  properties  of  these  high-order  methods,  combined  with  the  geometrical 
flexibility  afforded  by  the  use  of  DC  methods  on  unstructured  meshes,  can  lead  to  revolutionary 
advances  in  simulation  capability,  enabling  the  accurate  simulation  of  complex  phenomena  with 
a  wide  range  of  scales  from  first  principles,  while  reducing  discretization  errors  to  acceptable  lev¬ 
els.  Furthermore,  adjoint-based  sensitivity  analysis  techniques  were  also  included  and  investigated 
in  order  to  demonstrate  the  possibility  of  performing  simulations  with  quantifiable  error  bounds, 
the  enhancement  of  the  reliability  of  h-p  adaptive  methods,  and  enabling  the  solution  of  design 
optimization  problems  using  high-order  methods. 

The  Phase  I  portion  of  this  project  has  sought  to  develop  and  demonstrate  the  key  technologies 
and  capabilities  which  will  be  required  to  construct  a  production  level  simulation  tool.  The  produc- 
tionalization  of  these  capabilities  in  the  Phase  2  work  will  be  facilitated  by  making  extensive  use 
of  the  three-dimensional  unstructured  mesh  simulation  software  infrastructure  already  developed 
by  Scientific  Simulations  LLC  for  second-order  accurate  finite-volume  methods.  The  outcome  of  a 
successful  combined  Phase  1  and  2  award  will  be  a  revolutionary  high-order  simulation  capability 
which  will  be  marketed  to  government  and  commercial  clients  as  a  replacement  tool  to  current 
second-order  accurate  simulation  codes. 


2.  WORK  PERFORMED 

The  original  Phase  1  STTR  proposal  concentrated  on  several  specific  key  techniques  which  are 
required  to  be  demonstrated  in  order  to  justify  the  feasibility  of  the  overall  approach  for  constructing 
a  production  level  solver  in  subsequent  Phase  2  work.  Seven  individual  tasks  were  outlined  in  the 
Phase  1  proposal.  All  tasks  have  been  performed  during  the  course  of  this  project  and  results 
obtained  for  each  specific  task  are  reported  in  the  following  section.  The  seven  tasks  consist  of: 

•  Task  1:  Demonstrate  and  document  performance  of  Interior  Penalty  method  for  the  simu¬ 
lation  of  diffusion  phenomena  and  viscous  flows. 

•  Task  2:  Investigate  and  compare  the  effectiveness  of  shock  capturing  using  DC  methods 
with  limiters  versus  high-order  sub-cell  shock  capturing  with  artificial  dissipation  terms. 
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•  Task  3:  Extend  high-order  geometric  conservation  law  for  DG  methods  to  high-order  implicit 
Runge-Kutta  time-stepping  schemes. 

•  Task  4:  Demonstrate  adjoint-based  shape  optimization  for  aerodynamic  problems  based  on 
high-order  DG  discretizations  with  curved  boundaries. 

•  Task  5:  Demonstrate  goal-oriented  h-p  adaptivity  using  adjoint-based  error  estimation  for 
aerodynamic  problems  using  DG  discretizations. 

•  Task  6:  Extend  current  DG  solver  to  multi- physics  formulation  and  demonstrate  by  solving 
electromagnetics  problem. 

•  Task  7:  Validate  and  benchmark  current  DG  discretization  approach  with  existing  mature 
second-order  accurate  finite-volume  unstructured  mesh  solver. 

Additional  work  was  also  performed  on  other  tasks  which  were  considered  critical  to  the  demonstra¬ 
tion  of  the  feasibility  of  this  approach  for  large  scale  three-dimensional  simulations.  This  additional 
work  included  devising  a  technique  for  creating  curved  mesh  elements  in  boundary  layer  regions, 
and  extending  the  Navier-Stokes  Interior  Penalty  (IP)  discretization  to  three  dimensions. 

3.  RESULTS 

3.1  Task  1:  Interior  Penalty  Method  for  Diffusion 

The  original  discontinuous  Galerkin  discretization  developed  for  inviscid  flows  in  previous  work 

[1]  has  been  extended  to  the  full  Navier-Stokes  equations  using  the  interior  penalty  (IP)  method 

[2]  with  an  explicit  expression  developed  for  the  penalty  parameter  which  removes  the  previously 
heuristic  nature  of  this  approach  [3] .  This  approach  is  appealing  since  it  constitutes  one  of  the  sim¬ 
plest  DG  approaches  to  diffusion  terms,  retains  a  nearest-neighbor  stencil,  and  provides  equivalent 
performance  compared  to  more  complex  discretizations.  This  is  illustrated  in  Eigure  1  (reproduced 
from  Reference  [4] )  where  the  accuracy  achieved  for  the  solution  of  a  Poisson  equation  using  the  IP 
method  is  compared  with  the  more  commonly  used  BR2  DG  discretization  [5],  showing  equivalent 
accuracy  for  discretization  orders  of  accuracy  ranging  from  p=l  to  p=4.  The  discretization  of  the 
Navier-Stokes  equations  using  the  IP  method  has  been  implemented  for  both  triangular  and  quadri¬ 
lateral  elements,  while  the  efficiency  of  the  h-p  multigrid  solver  previously  developed  for  inviscid 
flows  has  been  enhanced  by  adding  a  line-implicit  solver  in  boundary  layer  regions  and  using  the 
entire  line-solver  driven  h-p  multigrid  algorithm  as  a  preconditioner  for  a  GMRES  krylov  method. 
Eigure  2  illustrates  the  solution  of  laminar  viscous  flow  (Re  =  5000)  over  a  slotted  airfoil  con¬ 
figuration  on  a  hybrid  mesh  including  triangular  and  quadrilateral  elements,  using  a  fourth-order 
accurate  DG  discretization  and  making  use  of  the  IP  method  for  the  discretization  of  the  viscous 
terms.  Eigure  2(c)  shows  the  convergence  rate  achieved  for  this  test  case  using  the  baseline  h-p 
multigrid  solver  with  an  element-implicit  (Jacobi)  smoother  on  each  multigrid  level,  a  line-implicit 
smoother,  and  using  the  latter  approach  as  a  preconditioner  for  GMRES.  The  inclusion  of  a  line- 
implicit  smoother  is  seen  to  be  crucial  for  obtaining  efficient  convergence  rates  in  the  presence  of 
anisotropic  meshes,  while  the  addition  of  GMRES  provides  further  modest  speedup. 
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Figure  1:  Comparison  of  accuracy  achieved  with  IP  discretization  and  BR2  discretization  for  Poisson  equation 
problem  in  two-dimensions  for  various  discretization  orders  of  accuracy. 


Figure  2:  (a)  Computed  Mach  contours  for  viscous  flow  over  slotted  airfoil  configuration  at  Mach  =  0.2  and 
Reynolds  number  of  5000  using  fourth-order  (p=3)  discretization  on  mixed-element  unstructured  mesh,  (b) 
Details  of  separation  pattern  obtained  for  this  test  case,  (c)  Convergence  of  h-p  multigrid  for  this  test  case  with 
element-implicit  solver,  line-implicit  solver,  and  used  as  a  preconditioner  for  GMRES. 
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(a)  (b) 

Figure  3:  (a)  Simulation  of  transonic  flow  over  NACA  0012  airfoil  at  Mach=0.8  and  1.25  degrees  incidence 
using  artificial  dissipation  approach  illustrating  sub-cell  shock  capturing  resolution  with  fourth-order  (p=3)  DG 
discretization,  (b)  Simulation  of  hypersonic  inviscid  flow  over  cylinder  at  Mach  =  6  using  third-order  (p=2)  DG 
discretization  with  limiting  procedure  for  robustly  capturing  strong  shock  wave. 

3.2  Task  2:  Shock  Capturing 

Two  approaches  have  been  investigated  for  shock  capturing  with  high-order  DG  discretizations.  The 
first  approach  is  based  on  the  use  of  artificial  dissipation  terms,  discretized  using  the  IP  method, 
to  smooth  out  and  capture  shocks.  A  captured  shock  for  an  inviscid  transonic  NACA0012  test  case 
is  shown  in  Figure  3(a),  showing  good  sub-cell  resolution  of  the  shock  wave.  While  this  approach 
works  well  for  relatively  weak  transonic  shocks,  it  has  been  found  to  be  problematic  for  strong 
shocks  encountered  in  hypersonic  flows.  Therefore,  a  non-linear  limiting  or  filtering  technique 
has  been  developed  for  capturing  strong  shocks.  The  approach  is  based  on  filtering  techniques 
originally  derived  for  image  reconstruction  problems  [6].  The  idea  is  to  use  this  approach  as  a 
filter  on  oscillatory  solutions  produced  near  under-resolved  or  discontinuous  phenomena  to  restore 
monotone  profiles.  Figure  3(b)  illustrates  the  solution  of  Mach  6  flow  over  a  cylinder  calculated 
using  the  non-linear  filtering  approach,  showing  a  crisp  nearly  monotone  shock  structure.  Results 
similar  to  those  displayed  in  Figure  3(a)  have  also  been  obtained  for  transonic  flow  problems.  These 
two  approaches  will  be  further  investigated  and  the  most  promising  strategy  will  be  used  in  the 
three-dimensional  solver  in  the  Phase  2  project. 

3.3  Task  3:  Arbitrary-Lagrangian-Eulerian  Formulations 

An  arbitrary  Langrangian-Eulerian  formulation  for  DG  discretizations  suited  for  solutions  on  dy¬ 
namically  deforming  meshes  which  guarantees  discrete  conservation  as  embodied  in  the  Geometric 
Gonservation  Law  (GGL)  has  previously  been  derived  for  first  and  second-order  time  discretiza¬ 
tions  [7].  In  the  current  work,  this  has  been  extended  to  higher-order  implicit  Runge-Kutta  time 
discretization  schemes.  For  a  first-order  backwards  difference  scheme,  the  basic  idea  consists  of 
formulating  the  problem  as  a  space-time  DG  discretization  with  constant  basis  functions  in  the 
time  direction.  The  extension  to  higher-order  IRK  schemes  is  then  obtained  by  applying  this  for¬ 
mulation  at  each  individual  stage  in  the  IRK  scheme.  This  formulation  has  been  validated  by  first 
verifying  that  uniform  flow  remains  an  exact  solution  in  the  presence  of  prescribed  mesh  deforma¬ 
tions,  including  the  deformation  of  curved  high-order  elements.  Subsequently,  it  was  verified  that 
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(a) 


(b) 


Figure  4:  (a)  Simulation  of  vortex  convection  problem  on  dynamically  deforming  mesh  for  testing  accuracy 
of  ALE  formulation.  (b)Temporal  accuracy  observed  for  static  and  moving  mesh  cases  with  various  temporal 
discretizations  for  GCL  compliant  ALE  DG  scheme. 

the  order  properties  of  third-  and  fourth-order  implicit  Runge-Kutta  schemes  are  preserved  in  the 
presence  of  moving  curved-element  meshes  for  the  test  case  of  a  convecting  inviscid  vortex  in  a 
uniform  flow.  These  results  are  summarized  in  Figure  4,  where  the  temporal  error  as  a  function  of 
time  step  size  is  compared  for  static  and  dynamic  mesh  cases. 

3.4  Task  4  :  Adjoint-Based  Shape  Optimization 

The  Phase  1  work  has  also  investigated  the  use  of  a  DG  adjoint  capability  for  driving  shape 
optimization  problems.  In  order  to  enable  the  computation  of  shape  sensitivities,  the  adjoint 
problem  must  be  first  formulated  and  solved.  The  adjoint  operator  is  formulated  by  first  forming 
the  complete  Jacobian  of  the  discretized  flow  equations  and  then  taking  the  transpose  of  this  matrix. 
This  constitutes  the  discrete  adjoint  approach,  and  is  easily  implemented  since  the  Jacobian  is 
available  from  the  linear  multigrid  solver  as  a  set  of  diagonal  block  sub-matrices,  which  are  easily 
transposed,  and  a  set  of  off-diagonal  block  sub-matrices,  which  are  transposed  and  swapped  to 
form  the  complete  transpose  operator.  Solution  of  the  adjoint  problem  is  then  performed  using  the 
same  h-p  multigrid  scheme  as  employed  for  the  flow  problem,  albeit  in  linear  form.  The  computed 
adjoint  field  must  then  be  multiplied  by  the  sensitivity  of  the  residual  with  respect  to  changes  in  the 
geometric  element  mappings,  which  may  include  the  effect  of  changes  in  element  curvature.  These 
terms  have  been  derived  and  demonstrated  for  an  inviscid  subsonic  airfoil  shape  optimization 
problem,  as  depicted  in  Figure  5.  Using  an  objective  based  on  the  difference  in  prescribed  and 
calculated  airfoil  surface  pressure,  the  target  surface  pressure  distribution  and  airfoil  shape  were 
recovered  in  40  design  cycles,  using  a  simple  steepest  descent  optimization  strategy  with  145  design 
variables  (Hicks-Henne  bump  functions)  used  to  define  the  airfoil.  The  target  airfoil  was  defined 
as  a  perturbation  of  the  initial  airfoil  which  ensured  the  realizability  of  this  optimization  problem. 
Although  this  represents  a  relatively  simple  optimization  problem,  it  serves  to  demonstrate  the 
ability  to  perform  shape  optimization  with  DG  discretizations  in  the  presence  of  curved  mesh 
elements.  This  demonstration  will  serve  as  the  basis  for  the  incorporation  of  a  shape  optimization 
capability  in  the  three-dimensional  solver  in  Phase  2,  where  the  principal  task  will  consist  of  deriving 
the  linearizations  of  the  residual  with  respect  to  the  element  mappings. 
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(a)  (b) 

Figure  5:  (a)  Evolution  of  designed  airfoil  from  initial  shape  to  final  shape  compared  with  target  airfoil  for  shape 
optimization  using  high-order  DG  discretization  of  the  Euler  equations  on  mesh  with  curved  surface  elements, 
(b)  Convergence  of  optimization  problem  as  measured  by  decrease  in  objective  function  versus  number  of  design 
cycles  using  steepest  descent  optimization  procedure. 

3.5  Task  5  :  Goal-Oriented  h-p  Adaptivity 

Adjoint-based  error-estimation  techniques  have  been  the  subject  of  considerable  interest  recently 
for  both  spatial  and  temporal  error  estimation  [8-11].  The  UW  group  has  investigated  adjoint 
formulations  for  both  finite-volume  and  DG  discretizations  over  the  last  several  years  [10-14].  The 
use  of  adjoint  methods  for  driving  h-p  refinement  strategies  with  DG  discretizations  has  been 
pursued  in  the  Phase  1  project.  Although  the  adjoint  error-estimation  strategy  remains  essentially 
unchanged  from  previous  implementations,  the  key  ingredient  for  enabling  h-p  refinement  is  a 
smoothness  indicator,  which  is  used  to  make  the  choice  between  h  or  p  refinement  in  regions  of 
large  error.  We  have  implemented  two  smoothness  indicators,  one  based  on  the  relative  size  of 
the  highest  solution  modes  versus  the  lower  solution  modes  [15],  and  another  based  on  the  size  of 
the  discontinuous  jumps  between  neighboring  elements  [16].  While  both  methods  give  satisfactory 
results,  the  second  approach  has  been  found  to  be  more  reliable  for  flows  with  strong  shocks. 
Figure  6(a)  illustrates  the  solution  obtained  for  hypersonic  flow  over  a  circular  cylinder  using 
h-p  adaptivity  with  the  jump-based  smoothness  indicator  for  the  adjoint-estimated  error  in  the 
surface  integrated  temperature  on  the  cylinder.  In  this  case,  extensive  h-refinement  occurs  in  the 
shock  region  where  the  solution  is  non-smooth,  and  p  refinement  occurs  in  the  region  between  the 
shock  and  the  cylinder  surface.  Note  that  no  refinement  (h  or  p)  occurs  in  areas  lateral  to  the 
cylinder,  since  these  areas  have  little  influence  on  the  objective  of  interest  (surface  temperature). 
Figure  6(b)  illustrates  the  convergence  of  the  simulation  objective  with  each  successive  refinement 
step,  showing  the  increasingly  accurate  error  prediction  provided  by  the  adjoint  procedure  at  each 
successive  refinement  level. 

3.6  Task  6:  Multi-Physics  Formulations 

One  of  the  consequences  of  the  local  compact  stencil  of  Discontinuous  Galerkin  discretizations  is 
that  a  simple  data  structure  consisting  of  elements  and  faces  can  be  used  to  support  the  discretiza¬ 
tion  of  various  types  of  equations.  In  the  Phase  1  work  we  have  demonstrated  this  by  extending  our 
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(a)  (b) 

Figure  6:  (a)  h-p  adaptive  refinement  pattern  for  hypersonic  flow  (Mach  6)  over  cylinder  based  on  error  in 
surface  integrated  temperature;(b)  Convergence  of  integrated  temperature  functional  with  refinement  levels  and 
error  prediction  provided  by  adjoint. 

existing  two-dimensional  DG  Euler  equation  solver  to  solve  the  governing  equations  for  electromag¬ 
netics.  Figure  7(a)  illustrates  an  inviscid  flow  solution  obtained  over  an  idealized  four-element  airfoil 
geometry,  while  Figure  7(b)  illustrates  the  solution  of  electromagnetic  scattering  over  a  cylinder 
using  the  same  code  modified  to  discretize  the  linear  Maxwell’s  equations.  A  further  extension  of 
this  code  to  include  linear  aeroacoustic  problems  is  currently  under  development.  In  this  case,  the 
discretization  is  closely  related  to  the  electromagnetic  problem,  with  the  exception  that  4  instead  of 
3  field  variables  are  required,  and  the  linear  equation  coefficients  depend  on  the  background  mean 
flow  field  which  must  be  read  in  at  run  time. 

3.7  Task  7:  Efficiency  Benchmarks 

The  performance  of  the  two-dimensional  Discontinuous  Galerkin  Navier-Stokes  solver  has  been 
compared  to  that  of  a  production  level  two-dimensional  unstructured  mesh  finite  volume  code 
(NSU2D  [17]).  For  these  comparisons,  the  line-implicit  h-p  multigrid  GMRES  solver  was  used  for 
the  high-order  Discontinuous  Galerkin  code,  using  discretization  orders  ranging  from  first-order  to 
fourth-order,  while  the  NSU2D  code  relies  on  a  second-order  accurate  finite- volume  vertex-based 
discretization,  and  makes  use  of  an  agglomeration  multigrid  technique  for  convergence  acceleration. 
The  two  codes  were  found  to  require  roughly  equivalent  computational  resources  for  cases  involving 
the  same  number  of  degrees  of  freedom,  thus  favoring  the  DG  scheme  due  to  the  higher  accuracy 
achieved  for  this  scheme  using  the  same  number  of  unknowns. 

3.8  Additional  Results 

An  additional  objective  addressed  during  this  project  relates  to  the  ability  to  generate  curved  mesh 
elements  which  conform  to  the  surface  geometry.  This  was  partly  motivated  by  comments  in  the 
debriefing  of  the  original  STTR  proposal  review,  but  also  became  necessary  upon  the  introduction 
of  stretched  meshes  near  the  surface.  In  our  approach,  additional  surface  points  within  each  surface 
element  are  obtained  by  accessing  the  original  spline  definition  of  the  geometry.  The  displacements 
between  these  new  points  and  their  projection  onto  the  linear  surface  edge  of  the  surface  elements 
is  then  used  to  define  the  curvature  of  the  surface  elements.  For  highly  stretched  meshes,  this  may 
lead  to  cross-over  with  neighboring  linear  elements,  and  therefore  the  curvature  is  propagated  to 
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(a) 
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Figure  7:  (a)  Inviscid  flow  solution  over  idealized  four  element  airfoil  configuration  using  p=4  (fifth-order)  accu¬ 
rate  DG  discretization.  Mach  =  0.2,  Incidence  =  0  degrees,  (b)  Instantaneous  x-component  magnetic  scattering 
field  computed  for  incident  radiation  over  circular  cylinder  using  p=4  (fifth-order)  accurate  DG  discretization. 


neighboring  elements  by  marching  along  lines  normal  to  the  boundary  layer  direction,  specifying 
the  same  displacements  on  each  subsequently  encountered  mesh  edge,  until  isotropic  elements  at 
the  edge  of  the  boundary  layer  are  encountered.  The  validity  of  the  resulting  curved  element  mesh 
is  then  confirmed  by  verifying  that  the  Jacobian  values  are  positive  at  all  quadrature  points  in  each 
element.  Note  that  it  is  not  sufficient  to  calculate  the  volume  of  these  elements,  as  this  corresponds 
to  an  average  Jacobian  value  at  all  quadrature  points.  Figure  8(a)  illustrates  the  curved  elements 
for  the  stretched  mesh  used  in  the  viscous  airfoil  flow  simulation  obtained  using  this  procedure. 
Extensions  to  three-dimensions  are  currently  underway,  using  the  VGRID  grid  generation  program 
with  a  facility  for  projecting  new  surface  points  onto  the  patched  geometry  description  used  by 
VGRID. 

The  previously  developed  three-dimensional  DG  Euler  solver  [18]  is  being  continuously  upgraded 
to  incorporate  the  advances  demonstrated  in  the  two-dimensional  setting  described  above.  Since 
this  will  ultimately  become  a  deliverable  production  code,  it  is  important  that  rigorous  software 
engineering  practices  be  adopted  from  the  outset  for  the  development  of  this  code.  To  this  end,  the 
three-dimensional  DG  solver  has  been  placed  under  the  SVN  software  version  control  system  [19] 
and  the  Hudson  [20]  regression  testing  system,  which  are  the  same  systems  used  in  the  HIARMS 
project  underway  in  the  University  of  Wyoming  research  group.  Individual  regression  tests  are  be¬ 
ing  implemented  and  added  to  the  test  suite  on  a  continual  basis.  The  planned  three-dimensional 
implementations  include  extension  to  the  Navier-Stokes  equations  using  the  IP  method,  support 
for  mixed  element  meshes,  and  the  generation  of  curved  surface  and  interior  elements  for  con¬ 
sistent  high-order  three-dimensional  meshes.  At  present,  a  three-dimensional  DG  Navier-Stokes 
discretization  using  the  analogous  IP  method  has  been  implemented  in  three-dimensions  for  tetra¬ 
hedral  meshes  with  straight  (non-curved)  elements  and  faces.  Figure  8(b)  illustrates  the  solution  of 
flow  over  a  backwards  facing  step  at  a  Reynolds  number  of  500  computed  with  the  element-implicit 
solver  in  the  absence  of  multigrid. 
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Figure  8:  (a)  Illustration  of  generation  of  curved  surface  and  interior  boundary  layer  elements  obtained  by 
propagating  curvature  information  from  boundary  along  lines  normal  to  boundary  layer,  (b)  Three-dimensional 
viscous  flow  computed  over  backwards  facing  step  using  fourth-order  (p=3)  DG  discretization  with  IP  method 
used  for  viscous  terms  on  fully  tetrahedral  mesh  for  subsonic  flow  at  Reynolds  number  of  500. 

4.  ESTIMATE  OF  TECHNICAL  FEASIBILITY 

The  Phase  1  project  has  been  designed  to  establish  the  technical  feasibility  of  constructing  a  robust 
high-order  accurate  Discontinuous  Galerkin  solver  for  complex  three-dimensional  configurations 
suitable  for  production  use.  The  specific  tasks  described  in  the  proposal  for  the  Phase  1  project 
were  chosen  to  demonstrate  specific  capabilities  which  will  be  required  for  the  construction  of  a 
successful  three-dimensional  production  simulation  capability. 

The  Phase  1  project  has  successfully  achieved  all  of  the  objectives  set  out  at  the  beginning  of 
the  project,  and  all  seven  technical  tasks  specified  in  the  proposed  work  plan  have  been  completed. 
The  Phase  1  work  has  demonstrated  the  feasibility  of  using  the  IP  method  for  discretizing  the 
Navier-Stokes  equations  in  two  and  three  dimensions.  Two  shock  capturing  schemes  have  been 
implemented  in  two  dimensions  and  used  successfully  to  capture  transonic  and  hypersonic  shocks. 
A  consistent  approach  for  formulating  high-order  DG  discretizations  in  the  presence  of  dynamically 
deforming  meshes  has  also  been  derived  and  demonstrated,  thus  making  feasible  the  application  of 
this  simulation  approach  to  more  complex  problems  with  moving  or  deforming  bodies.  An  adjoint 
solution  technique  has  been  developed  and  used  to  demonstrate  the  potential  for  shape  optimization 
in  the  presence  of  curved  surface  elements  with  high-order  accuracy,  as  well  as  for  predicting  error 
in  specific  functional  objectives  and  for  driving  h-p  adaptive  schemes,  which  hold  great  promise  for 
delivering  highly  accurate  and  certifiably  correct  solutions.  Using  the  same  software  framework  and 
data-structures,  the  two-dimensional  DG  code  was  extended  to  electromagnetic  problems,  illustrat¬ 
ing  the  potential  for  solving  multi-physics  problems  within  a  single  software  framework  based  on 
DG  discretizations.  A  technique  for  handling  stretched  meshes  in  the  presence  of  curved  surface  ge¬ 
ometries  has  also  been  devised,  and  strategies  for  extending  this  approach  to  three  dimensions  have 
been  developed.  Finally,  the  two-dimensional  DG  Navier-Stokes  solver  was  benchmarked  against 
a  production  finite- volume  Navier-Stokes  solver  and  found  to  be  competitive  in  terms  of  required 
computational  resources. 
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These  findings  validate  a  set  of  important  capabilities  which  will  be  instrumental  in  the  con¬ 
struction  of  a  three-dimensional  production  solver  based  on  DG  discretizations.  At  present,  a 
three-dimensional  Navier-Stokes  solver  based  on  high-order  DG  discretizations  and  using  the  ef¬ 
ficient  h-p  multigrid  solver  originally  developed  in  two  dimensions  has  been  developed  and  has 
been  validated  on  simple  geometries.  The  successful  demonstration  of  these  tasks  in  the  Phase  1 
project,  along  with  the  current  operational  status  of  our  three-dimensional  solver,  and  currently 
planned  extensions,  present  a  strong  case  for  the  feasibility  of  achieving  the  goals  set  out  in  the 
Phase  1  proposal,  namely  the  construction  of  a  high-order  accurate  simulation  tool,  suitable  for 
complex  geometries,  usable  in  the  presence  of  dynamically  deforming  geometries,  and  capable  of 
incorporating  the  latest  advances  in  numerical  simulation  technology,  such  as  adaptive  refinement 
and  design  optimization. 


5.  CONCLUSION  AND  FURTHER  WORK 

The  tasks  completed  in  the  Phase  1  project  were  mostly  carried  out  in  the  two-dimensional  setting. 
Work  is  currently  underway  to  extend  these  capabilities  into  the  three-dimensional  Navier-Stokes 
DG  code  which  has  been  developed  in  the  project  and  described  above.  In  order  to  incorporate 
all  the  capabilities  demonstrated  in  the  Phase  1  project,  a  re-engineering  of  the  original  three- 
dimensional  code  is  required.  The  three-dimensional  code  must  be  capable  of  solving  multiphysics 
problems,  handling  cases  with  dynamically  deforming  meshes,  where  individual  mesh  elements 
may  have  different  orders  of  geometric  mappings  (i.e.  straight-faced  elements  employing  linear 
mappings,  curved  elements  employing  higher-order  mappings).  Individual  mesh  elements  may  also 
have  different  solution  accuracies  (p  orders)  in  order  to  enable  h-p  adaptive  refinement  schemes. 
Thus,  a  new  data-structure  has  been  devised  where  mesh  elements  may  have  different  shapes  (mixed 
element  meshes),  may  have  individually  specifiable  solution  order,  individually  specifiable  geometric 
mapping  orders,  and  where  all  ALE  terms  are  incorporated  and  high  enough  quadrature  rules  are 
implemented  in  order  to  guarantee  conservation  in  the  presence  of  moving  meshes.  This  data 
structure  is  currently  being  implemented  in  a  linked  list  fashion,  and  will  enable  all  the  capabilities 
demonstrated  in  the  Phase  1  project  to  be  integrated  into  the  three-dimensional  solver. 

Based  on  the  demonstrated  results  from  the  Phase  1  project,  it  is  anticipated  that  continued 
development  of  this  technology  and  extension  into  the  three-dimensional  setting  following  the  path 
described  above,  will  revolutionize  our  ability  to  simulate  complex  engineering  problems  with  high 
fidelity  in  the  near  future. 
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11 


FA9550-09-C-0021:  High-Order  Modeling  of  Applied  Multi-Physics  Phenomena 

Scientific  Simulations  LLC 


References 

[1]  C.  Nastase  and  D.  J.  Mavriplis.  High-order  discontinuous  Galerkin  methods  using  a  spectral  multigrid 
approach.  Journal  of  Computational  Physics,  213(l):330-357,  March  2006. 

[2]  D.  N.  Arnold.  An  interior  penalty  finite  element  method  with  discontinuous  elements.  SIAM  Journal 
of  Numerical  Analysis,  19:742-760,  1982. 

[3]  K.  Shahbazi.  An  explicit  expression  for  the  penalty  parameter  of  the  interior  penalty  method.  Journal 
of  Computational  Physics,  205(2) :40 1-407,  May  2005. 

[4]  D.  J.  Mavriplis,  C.  Nastase,  L.  Wang,  K.  Shahbazi,  and  N.  Burgess.  Progress  in  higher-order  discon¬ 
tinuous  Galerkin  methods  for  aerospace  applications.  AIAA  Paper  2009-0601,  presented  at  the  47th 
Aerospace  Sciences  Meeting,  Orlando,  FL,  January  2009. 

[5]  F.  Bassi  and  S.  Rebay.  GMRES  discontinuous  Galerkin  solution  of  the  compressible  Navier-Stokes 
equations.  Discontinuous  Galerkin  Methods:  Theory,  Gomputations  and  Applications,  K.  Gockburn, 
W.  Shu  (Eds.),  Springer,  Berlin,  2000. 

[6]  L.  I  Rudin,  S.  Osher,  and  E.  Fatemi.  Nonlinear  total  variation  based  noise  removal  algorithms.  Physica 
D,  60:259-268,  1992. 

[7]  D.  J.  Mavriplis  and  G.  Nastase.  On  the  geometric  conservation  law  for  high-order  discontinuous  Galerkin 
discretizations  on  dynamically  deforming  meshes.  AIAA  Paper  2008-0778,  January  2008. 

[8]  J.  T.  Oden  and  S.  Prudhomme.  On  goal-oriented  error  estimation  for  elliptic  problems:  Application  to 
control  of  pointwise  errors.  Comput.  Methods  Appl.  Mech.  Engr.,  176:313-331,  1999. 

[9]  D.A.  Venditti  and  D.L.  Darmofal.  Grid  adaptation  for  functional  outputs:  Application  to  2-D  inviscid 
compressible  flow.  Journal  of  Computational  Physics,  176:40-69,  2002. 

[10]  L.  Wang  and  D.  J.  Mavriplis.  Adjoint-based  h-p  adaptive  discontinuous  Galerkin  methods  for  the 
compressible  Euler  equations.  AIAA  Paper  2009-0952,  January  2009. 

[11]  K.  Mani  and  D.  J.  Mavriplis.  Error  estimation  and  adaptation  for  functional  outputs  in  time-  dependent 
fiow  problems.  AIAA  Paper  2009-1495,  presented  at  the  AIAA  Aerospace  Sciences  Meeting,  Orlando 
FL,  January  2009. 

[12]  D.  J.  Mavriplis.  Discrete  adjoint-based  approach  for  optimization  problems  on  three  dimensional  un¬ 
structured  meshes.  AIAA  Journal,  45(4):741-750,  April  2007. 

[13]  K.  Mani  and  D.  J.  Mavriplis.  Unsteady  discrete  adjoint  formulation  for  two-dimensional  fiow  problems 
with  deforming  meshes.  AIAA  Journal,  46(6):1351-1364,  2008. 

[14]  D.  J.  Mavriplis.  Solution  of  the  unsteady  discrete  adjoint  for  three-  dimensional  problems  on  dynamically 
deforming  unstructured  meshes.  AIAA  Paper  2008-0727,  January  2008. 

[15]  P.  Persson  and  J.  Peraire.  Sub-cell  shock  capturing  for  discontinuous  Galerkin  methods.  AIAA  Paper 
2006-0112,  January  2006. 

[16]  L.  Krivodonova,  J.  Xin,  J.  F.  Remade,  N.  Ghevaugeon,  and  J.  Flaherty.  Shock  detection  and  limiting 
with  discontinuous  Galerkin  methods  for  hyperbolic  conservation  laws.  Applied  Numerical  Mathematics, 
48:323-338,  2004. 

[17]  W.  O.  Valarezo  and  D.  J.  Mavriplis.  Navier-Stokes  applications  to  high-lift  airfoil  analysis.  AIAA  J.  of 
Aircraft,  23(3):457-688,  1995. 


12 


FA9550-09-C-0021:  High-Order  Modeling  of  Applied  Multi-Physics  Phenomena 

Scientific  Simulations  LLC 


[18]  C.  Nastase  and  D.  J.  Mavriplis.  A  parallel  hp-multigrid  solver  for  three-dimensional  discontinuous 
Galerkin  discretizations  of  the  Euler  equations.  AIAA  Paper  2007-0512,  January  2007. 

[19]  Subervsion  website,  http://subversion.tigris.org/. 

[20]  Hudson  website,  https://hudson.dev.java.net/. 


13 


